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Abstract 

We study the uncertainties in parton distributions, determined in global fits to deep 
inelastic and related hard scattering data, due to so-called theoretical errors. Amongst 
these, we include potential errors due to the change of perturbative order (NLO NNLO), 
ln(l/x) and ln(l — x) effects, absorptive corrections and higher-twist contributions. We 
investigate these uncertainties both by including explicit corrections to our standard global 
analysis and by examining the sensitivity to changes of the x, Q^, W"^ cuts on the data that 
are fitted. In this way we expose those kinematic regions where the conventional DGLAP 
description is inadequate. As a consequence we obtain a set of NLO, and of NNLO, 
conservative partons where the data are fully consistent with DGLAP evolution, but over 
a restricted kinematic domain. We also examine the potential effects of such issues as 
the choice of input parameterization, heavy target corrections, assumptions about the 
strange quark sea and isospin violation. Hence we are able to compare the theoretical 
errors with those uncertainties due to errors on the experimental measurements, which 
we studied previously. We use W and Higgs boson production at the Tevatron and the 
LHC as explicit examples of the uncertainties arising from parton distributions. For many 
observables the theoretical error is dominant, but for the cross section for W production 
at the Tevatron both the theoretical and experimental uncertainties are small, and hence 
the NNLO prediction may serve as a valuable luminosity monitor. 



^ Royal Society University Research Fellow. 



1 Introduction 



A knowledge of the partonic structure of the proton is an essential ingredient in the analysis of 
hard scattering data from pp or pp or ep high energy collisions. However, only the (or scale) 
dependence of the parton distributions can be calculated from perturbative QCD. Perturbative 
QCD cannot fix their absolute values. Moreover, non- perturbative techniques arc still far from 
being able to predict reliable magnitudes. Rather, to determine the distributions, it is necessary 
to resort to global analyses of a wide range of deep inelastic and related hard scattering data. 
The Bjorken x dependence of the distributions is parameterized at some low scale, and a fixed 
order (either LO, NLO or NNLO) DGLAP evolution performed to specify the distributions at 
the higher scales where data exist. Much attention is now being devoted to obtaining rehable 
uncertainties on the parton distributions obtained in this way. One obvious uncertainty is due 
to the systematic and statistical errors of the data used in the global fit. We will call these the 
experimental errors on the parton distributions and on the physical observables predicted from 
them. In fact, these so-called experimental uncertainties of the partons have so far been the 
main focus of attention. They have been estimated by several groups [1]-[10], working within a 
NLO framework using a variety of different procedures. For instance, in a previous paper [10] 
we estimated the parton errors using a Hessian approach in which we diagonalised the error 
matrix and then used the linear propagation of errors to estimate the uncertainty on a variety 
of typical observables. We confirmed that this simple approach works well in practice by using 
a more rigorous Lagrange multiplier method to determine the errors on the physical quantities 
directly. We also compared our results with those obtained in similar analyses performed by 
the CTEQ collaboration [4, 5, 6]. 

Besides the experimental errors, there are many other sources of imcertainty associated with 
the global parton analysis. These are the concern of the present paper. They may loosely be 
called theoretical errors. That is, we use theoretical errors to denote all uncertainties on the 
predicted observables other than those that arise from the systematic and statistical errors of 
the data that are included in the global fit. The various theoretical errors may be divided 
into four categories. Uncertainties due to (i) the selection of data fitted, (ii) the truncation of 
the DGLAP perturbation expansion, (in) specific theoretical effects^ (namely In ln(l — x), 
absorptive and higher twist corrections) and (iv) input assumptions (such as the choice of 
parameterization, heavy target corrections, whether isospin violation is allowed and the form 
of the strange quark sea). These theoretical errors are discussed in turn in the following four 
sections. Then, in Section 6, we study the implications of the error analysis for determining the 
value of as(M|) from the global fit, and for the accurate predictions of particularly relevant 
observables at the Tevatron and the LHC. 

^Of course, differences can also arise from alternative methods of treating the behaviour near the heavy flavour 
thresholds. However, this issue is now well understood. Either fixed-flavour-number schemes or, preferably, one 

of the; sdcction of variable- flavour- number schemes can be used for a correct description [11]. Different choices 
will result in slight differences in the extracted partons, but only minimal variation in predictions for physical 
quantities. 
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2 Uncertainties due to the selection of data fitted 



In principle, if the DGLAP formalism is valid and the various data sets are compatible, then 
changing the data that arc included in the global analysis should not move the predictions 
outside the error bands. In practice this is not the case. As an extreme example, consider the 
analysis of the HI and BCDMS deep inelastic data for F2 [7]. The absence of the other fixed 
target data and, in particular, the Tevatron jet data from the global parton fit results in a 
gluon which is much smaller for x > 0.3, and hence larger at small x, from those obtained in 
global fits to a wider range of data. Another example occurs in the analysis of Ref . [2] , where 
fitting to a small subset of the total data (namely Hl(94) [12], BCDMS [13] and E665 [14] data 
for F2) yields a surprisingly low value as — 0.112 ± 0.001, and a very hard high x down-quark 
distribution. Even fits of different groups, at the same perturbative order, to basically the same 
data lead to unexpectedly sizable differences in partons and in predictions for observables. For 
example, compare the MRST2002 and CTEQ6 predictions at NLO for the cross sections for W 
and Higgs boson production at hadron colliders, shown in Fig. 15 of Ref. [10]. We see that the 
CTEQ6 values at the Tevatron are more than 1.5% smaller for W production and 8% smaller 
for Higgs production; differences which are larger than expected from an analysis based on the 
experimental errors of the data used in the fits. 

We begin by investigating the stability of the global parton analysis to different choices 
of the data cuts (VFcut, ^^cut, Qcut)) defined such that data with values olW^xoxQ^ below 
the cut are excluded from the global fit, with the implicit assumption that the remaining data 
can be described by pure leading-twist DGLAP evolution. We find the minimum values of the 
data cuts for which this stability occurs. The kinematic variable W is the invariant mass of 
the system X recoiling against the scattered lepton in deep inelastic lepton-proton scattering 
Ip ^ IX. It follows that 

W"^ c^Q'^{l-x)/x. (1) 

In the remainder of this section, all fits are performed at next-to-leading order (NLO) unless 
otherwise stated. In Section 3 we discuss what happens when we use the NNLO formalism. 

2.1 Effect of the cut on W'^ 

In the original MRST global analyses we have fitted to data with W"^ > 10 GeV^, assuming 
that this was sufficiently high to avoid resonance structure, large ln(l — x) effects and associated 
higher-twist corrections. However, we had no quantitative justification for this precise choice, 
and, noting that it resulted in a systematically poor fit to SLAG data with W'^ > 10 GeV^, we 
subsequently raised the cut to W^^^ — 12.5 GeV^ [15]. This provided an acceptable description 
of the SLAG data. In order to make a more systematic investigation of the stability of the fit 
to varying this cut, we performed a series of global analyses with VF^ut ranging from 12.5 to 
25 GeV^. When raising the cut from 12.5 to 15 GeV^ we find to the remaining data improved 
by only 4, while an increase of W^^^ from 15 to 20 or 25 GeV^ resulted in no significant 
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improvement. Taking these results at face value, we conclude that W^^^. — 15 GeV^ is a 
conservative choice and that there is no reason to eliminate even more data. 

However, inspection of the description of the SLAG and BCDMS data in this low W'^ 
region shows a lack of compatibility of the two data sets in the region where they overlap; 
specifically for in the interval (6, 15) GcV^ for the higher values of x. Thus the stability 
achieved at W^cut ~ GeV^ corresponds to the SLAG data completely disappearing from the 
fit. Hence, while the resulting fit describes the BGDMS data well, an extrapolation to only 
slightly lower values gives a very poor description of the SLAG data. This implies that the 
achieved stability is an artifact of the incompatibility of the two data sets for Indeed, 
when phenomenological higher-twist contributions are introduced into the analysis they still 
have significant impact for W'^ > 15 GeV^, see Fig. 2 of Ref. [16]. This implies that a genuinely 
conservative choice of cut would be Wcut ~ 20 GeV^. However, because a good description is 
obtained of the only available data in the (15, 20) GeV^ interval we may set Wcut = 15 GeV^ 
without prejudicing the analyses. Future measurements of Fj* at HERA in this kinematic 
domain would clearly be valuable. 

2.2 Effect of the choice of the cut on x 

In Table 1 we show the values of for global analyses performed for different values of Xcut, 
together with the number of data points fitted. Each column represents the values corre- 
sponding to a fit performed with a different choice of the cut in x. For example, if only data 
above Xcut — 0.001 are fitted, then the total — 2119, with a contribution of — 2055 coming 
from the subset of data with x > 0.0025, and — 2012 from the subset with x > 0.005, etc. 

To obtain a measure of the stability of the analysis to changes in the choice of Xcut, we 
compare fits in adjacent columns, that is with {xcut)i+i and (a;cut)i- In particular, it is infor- 
mative to compare the contributions to their respective values from the subset of data with 
X > {xcut)i+i- If stability were achieved, then we would expect the difference Ax^ between 
these two contributions to be very small. We stress that these two x^ contributions describe 
the quality of the two fits to the same subset of the data. Thus, as we shall explain below, 
a measure, A-"*"^, of the stability of the analysis is A^^ divided by the number of data points 
omitted when going from the fit with (.Xcut)i to the fit with {xcut)i+i- For example, if we raise 
the Xont from 0.001 to 0.0025 then Ax^ = 2055 - 2040 for the data with x > 0.0025, and the 
number of data points omitted is 1961 - 1898 = 63. Thus the measure A[]:[][]f = 15/63 = 0.24, 
as shown in the last row of Table 1. If we were to start from a fit with Xcut below the value at 
which the theoretical framework is expected to be valid (due to neglected Inl/x effects, etc.), 
then we would expect A-"*"^ to be significantly non-zero. As Xcut is increased, then A^"^^ should 
decrease and, in the ideal case, approach and remain near zero; indicating that we are in the 
stable domain where the theoretical framework is appropriate. In this way the behaviour of 
A*"*"^, as Xcut is changed, acts as an indicator of the stability of the analysis. 
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X^{x > 0.01) 


1965 


1961 


1953 


1934 


1917 


1916 




0.19 0.10 0.24 0.28 0.02 



Table 1: Each column shows the values obtained from a NLO global analysis with a different 
choice of Xcut, together with the number of data points fitted and the value of as{Mz) obtained. 
The first entry in a given column is the total x^, and the subsequent entries show the 
contributions to x^ from subsets of the data that are fitted. The quantity A-^^, shown in the 
final row, is a measure of stability to changing the choice of Xcut, as explained in the text. In 
these analyses we take the default cut in Q^, that is Qcut — 2 GeV^. 

Inspection of the values of A^"*"^ in the last row of Table 1 shows a significant improvement 
in the quality of the fit each time Xcut is raised by an amount corresponding to the omission of 
about a further 70 data points, until the final step when Xcut is increased from 0.005 to 0.01, 
when we see that there is no further improvement at all. In fact, raising Xcnt from 0.01 to 0.02 
confirms this stability. Indeed, A*^^ is greater when raising x^ut from 0.0025 to 0.005 than in 
any of the previous steps. Hence we conclude that x ~ 0.005 is a safe choice of Xcnt- Below this 
value there is a degree of incompatibility between the data and the theoretical description. At 
each step the gluon distribution extrapolated below Xcut becomes increasingly smaller, allowing 
the higher x gluon to increase (and to carry more momentum), see Fig. 1. Because almost no 
momentum is carried by the gluon at very small x, this redistribution of the gluon momentum 
becomes increasingly possible as Xcut is raised, hence explaining why A^"*"^ has a tendency to 
increase, until stability is finally reached. In general, we see that there is a slight decrease of 
asiM^r) from our standard value of 0.1197 to 0.1178. 

2.3 Effect of the choice of the cut on 

Repeating the above procedure, but now performing fits with different choices of the cut, we 
obtain the results shown in Table 2. Here the behaviour of the values of A*'^^ is not so dramatic 
as for the study of x^ut- As a consequence it is more difficult to select the value of Ql^^ for which 
the theoretical description becomes safe. However, there is a general decrease in the value of 
A*"^^ as the value Qcut i^ increased. Certainly the choice Qcut < GeV^ appears inappropriate. 



4 




Figure 1: The gluon distribution obtained in NLO global fits with different values of Xcut, that 
is Xcut taken to be 0.0002 (dashed curve), 0.001 (dotted) and 0.005 (dot-dashed), compared to 
the default set with Xcut — (continuous curve). 
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Table 2: Each column shows the values obtained from a NLO global analysis with a different 
choice of the cut in Q^, together with the number of data points fitted and the value of as{Mz) 
obtained. The first entry in a given column is the total x^, and the subsequent entries show 
the contributions to x^ from subsets of the data that are fitted. A*"*"^, shown in the final row, 
is a measure of stability to changing the choice of Qcut; ^ explained in the text. 

and Qcut > 14 GeV^ is definitely acceptable. We therefore take the reasonably conservative 
choice of Ql^^ — 10 GeV^. This gradual reduction of A^"*"^ is an indication of the presence 
of higher-order corrections, whose relative strength falls off only hke l/lnQ^. If higher-twist 
corrections were the dominant effect, then we would expect a more dramatic reduction of A*"^^ 
at some low value of Qcut- 

2.4 Effect of a cut on the product xQ"^ 

The theoretical uncertainties at small x and small may be strongly correlated. That is, the 
main theoretical uncertainty at small x is due to higher-twist effects rather than higher-order 
contributions. In order to investigate this possibility, we perform a series of global fits, each 
with a different choice of cut on the product xQ"^ . We begin with (.T(5^)c„t = 0.001 GeV^, which 
corresponds to a loss of 42 data points, and proceed in steps of a loss of about a further 50-100 
data points until stability is achieved. This occurs when (a;(5^)cut — 0.6 GeV^, at which point 
589 data points have been removed. The corresponding value of q;5(M|) = 0.1183, and the 
gluon distribution, are both very similar to those obtained with Xcvx — 0.005. Hence we arrive 
at a similar final result to that of applying an Xcnt alone, but with the loss of about twice as 
much data. Indeed we have removed all HERA data below x ~ 0.003, and have only a handful 
of high points left for x < 0.005. Thus to cut on xQ^ appears to be much more inefficient 
than to cut on x alone. 
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2.5 Combined cut on x and 

Finally we check the effect of imposing a range of cuts on x having already chosen the con- 
servative cut of Qcut = 10 GeV^, and alternatively choosing different cuts in having 
already selected the safe Xcut of 0.005. In each case the additional cut leads to further improve- 
ments in the quality of the fit, and complete stability is only achieved when the combined cuts 
a^cut = 0.005 and = 10 GeV^ are imposed. (A*"*"^ = 0.13 when going from Xcut = 0.0025 
to 0.005 at Ql^^ = 10 GeV^ and A^+^ = 0.05 when going from Ql^, = 7 GeV^ to 10 GeV^ 
at Xcut = 0.005). This combined cut yields, in our opinion, parton distributions largely free 
from theoretical uncertainties, but only within this restricted kinematic domain. We denote 
this conservative parton set by MRST(cons). The corresponding value of as{M^) is 0.1162, 
and the gluon distribution is similar in shape to that obtained for the fit with a;cut = 0.005 
(and Qcut = 2 GeV^). To estimate the error of this MRST(cons) prediction of ag we use a 
tolerance of Ax^ = 5, instead of our previous choice of Ax^ = 20 [10], since we now have a more 
compatible collection of data, without so many accompanying 'tensions' between different data 
sets. However, the reduced amount of fitted data with, in particular, a more limited range of 
, yields an intrinsically larger error so that our final uncertainty on as is again approximately 
±0.002 (despite the smaller tolerance). 

Fig. 2 compares the MRST(cons) partons distributions with those of our default set. We 
stress that these partons have a limited range of applicability. In order to investigate the 
theoretical uncertainty outside the kinematic domain, x > 0.005 and Q"^ > 10 GeV^, of the 
fitted data, we next study a variety of possible theoretical corrections to the NLO leading-twist 
framework. 

3 Uncertainty associated with change from NLO to NNLO 
fit 

Clearly there will be uncertainties associated with the truncation of the DGLAP evolution at a 
given perturbative order. In the past, these have usually been estimated by seeing the effect of 
changes of scale on physical observables in the NLO analysis. However, this approach is flawed 
since it provides no information about higher logs in 1/x and (l—x) at higher orders, and now 
we can do better. The deep inelastic coefficient functions are known at NNLO [17]. Moreover, 
valuable partial information has been obtained about the NNLO splitting functions [18, 19]. 
This greatly limits the possible behaviour of the splitting functions down to quite small values 
of X. Indeed, van Neerven and Vogt [20] have constructed a range of compact analytic functions 
that are all compatible with the available information. 

3.1 Effect of NNLO corrections 

We have performed NNLO global analyses in two previous publications [21, 22]. The dominant 
effects of the NNLO corrections are indeed the additional ln(l — x) terms in the nonsinglet 
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Figure 2: Comparison of MRST(cons) partons with the default NLO set MRST(2002). The 
former partons are only reliable for x > 0.005 and > 10 GeV^. 
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coefficient functions wfiicfi influence large x, and the behaviour of the coefficient functions 
and splitting functions at small which is heavily determined by the leading \u{l/x) terms. 
The additional terms in the nonsinglet coefficient function lead to an enhancement of the 
large x structure functions, and hence a small, but significant, reduction in the valence quark 
distributions at large x. However, because this enhancement is proportional to q;|((5^) it falls 
quickly with and the evolution of the structure functions increases at large x. This enables 
a slightly better fit to be obtained, but the natural increase in evolution speed results in a 
lower value of ^^(Ml), i.e. it falls from ~ 0.119 at NLO to ~ 0.116 at NNLO, which is mainly 
determined by the high-x structure function evolution. 

At small X the speed of evolution of -^2(0;, Q'^) is also increased, both by the behaviour of the 
NNLO splitting function P^^\x) and the NNLO coefficient function C2g{x). The evolution of 
the gluon is decreased, however, due to the leading \n{l/x)/x term in P^^J{x) having a negative 
coefficient. In our first analysis we found that the small x gluon at NNLO was far below that at 
NLO, but this was largely due to the initial estimate of P^{x), which was reduced significantly 
when some additional moments became available. In our more recent analysis we do indeed 
find that the NNLO gluon is a little smaller at small x due to the natural increase in evolution 
of F2{x,Q'^). However, the large x gluon needs to be larger at NNLO because the decrease in 
the large x quarks discussed above reduces the high-E'r Tevatron jet cross section, which must 
be corrected by an increase in the large-x gluon distribution.^ This redistribution of the gluon 
between large and small x is qualitatively in agreement with the momentum sum rule, which 
is a strong constraint on the form of the gluon. Indeed, the complete NNLO fit works nicely, 
and displays a very slight improvement in quality compared to NLO. 

However, it is not clear whether additional large logarithms in {1 — x) and 1/x beyond 
NNLO will lead to further modifications to the partons, and any resulting predictions. We 
have produced predictions for W and Z production at the Tevatron and the LHC in [21, 22], 
and also for the longitudinal structure function Fi{x,Q'^) (the NNLO coefficient function for 
this having been estimated [21] in the same way as the NNLO splitting functions). The former 
are quite stable, showing changes of ~ 4% when going from NLO to NNLO. However, this is 
larger than the ~ 2% 'experimental' error estimated at NLO [6, 10], and is dependent on the 
quark distributions, which are directly obtained by comparison with the structure functions. 
In contrast, Fl{x, Q^) is dependent mainly on the gluon which has no direct constraint at low 
X, so that when one goes from LO to NLO to NNLO the variation is very large, i.e., ~ 20% at 
X = 0.0001, even for > 100 GeV^. This variation is driven largely by the ln(l/a;) terms, and 
it is not clear if any sort of convergence has been reached at NNLO. The same is potentially 
true for other gluon dominated quantities. 

^There is not, at present, a full NNLO calculation of the jet cross section available, but a calculation of 
leading threshold logarithms [23] suggests that the NNLO contribution is not large or very E'r-dependent. 
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3.2 Cuts on Data at NNLO 



We might hope that at NNLO the cuts on data required to achieve stabihty are rather less 
stringent than at NLO. However, in the same way that the quahty of the global fit only improves 
very slightly going from NLO to NNLO, the size of the cuts in W'^, and x does not change 
much. Nevertheless, there are some advantages to be gained at NNLO as we discuss in detail 
below. 

Raising W^^^ from 12.5 GeV^ has a similar effect as at NLO. However, this stability at 
VFj,^^. = 12.5 GeV^ or at most 15 GeV^ is subject to the same caveat as at NLO, i.e. a certain 
incompatibility of data for ~ 15 GeV^. It is also true that the low W"^ data is one of 
the major constraints on as{M],), which is rather lower at NNLO than NLO, being more 
consistent between different data sets at NNLO. Also, if one extrapolates the NNLO fit below 
— 12.5 GeV^ the departure between data and theory occurs more slowly than at NLO. 
Hence, although the 'optimal' value of FF^fut does not change at NNLO there are indications 
that the theoretical description is improving at low W"^. We will return to this point later. 

At NNLO there is a significant improvement each time x^ut is raised until x^ut = 0.005, as at 
NLO. However, the improvement in going from x^ut = 0.001 to .Tcut = 0.005 is not as great as at 
NLO. The dominant improvement in the quality of the fit is due to a much improved fit to the 
Tevatron jet data due to a readjustment of the high x gluon. This was also the case at NLO, 
but there was a far more significant improvement in the description of the NMC data and the 
X > 0.005 HERA structure function data in that case. As at NLO, each time the value of Xcut is 
raised there is a reduction in the very low x gluon, and a corresponding general increase in the 
high and moderate x gluon. However, the effect is far less pronounced at NNLO than it was at 
NLO, as seen in Fig. 3 which compares the default gluons to the Xcut — 0.005 versions at both 
NLO and NNLO. Similarly, while q;5(M|) decreased as Xcut was lowered at NLO, it is almost 
completely insensitive to a^cut at NNLO. Hence, although the value of x above which we have 
complete confidence in our partons is the same at NNLO as at NLO, the changes compared to 
the default set, and hence the uncertainties involved with the default set, are much reduced at 
NNLO. 

It is a similar story for cuts. As at NLO there is a continual improvement until Qcut — 
10 GeV^, above which stability sets in, and as at NLO this improvement is relatively gradual 
between Q^ut — 2 GeV^ and Qcut — 10 GeV^, suggesting that higher twist is not most obviously 
the remaining correction to theory required. However, as with x cuts, the total improvement 
in the quality of the fit and the degree of readjustment of the partons is not quite as large as 
at NLO. 

Finally, we investigate the effects of cuts on both x and Q^. With our experience at NLO 
we would expect stability to be achieved when both Xcut and Qcut ^'^^^ their respective 
individual values for stability. However, there is a slight improvement at NNLO as compared to 
NLO. In this case A*"*"^ is negligible when going from Ql^^ = 7 GeV^ to 10 GcV^ at Xcut = 0.005. 
However, A-"*"^ = 0.12 when going from Xcut = 0.0025 to 0.005 at Ql^^ = 7 GeV^, and similarly 
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Figure 3: Comparison of MRST gluons obtained from Xcut = 0.005 with the default gluons 
MRST(2002) at both NLO and NNLO. 
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if Ql^^ — 10 GeV^. Hence, our conservative set of NNLO partons is obtained with the shghtly 

less severe Q^^^ = 7 GcV^ and x^ut = 0.005. We also note that the modification of the 
gluon distribution for the MRST(cons) set is far smaller at NNLO than at NLO, being similar 
to the Xcnt = 0.005 alone set as shown in Fig. 3, and the value of as{M^) reduces only to 
0.1153 ± 0.002, a much smaller reduction than at NLO. As at NLO, the error corresponds to 
the reduced tolerance Ax^ = 5. The detailed comparison of the 'conservative' set with the 
default set of partons is shown in Fig. 4, and one sees that the deviation is rather smaller than 
at NLO, particularly at small x. 

Hence, the investigation of the full range of cuts on the data that are fitted results in a set of 
'conservative' NNLO partons which are only completely reliable over a slightly extended range 
compared to NLO. However, the change in these partons compared to the default set within 
this range, and more especially outside the range, is much smaller than at NLO (and the same 
is true for the change in as{M'^)). This implies that predictions made using NNLO partons are 
considerably more reliable than those using NLO partons, even before we consider the extra 
precision obtained simply by using the expressions for cross sections at a higher order. We will 
discuss this in more detail later. 

4 Specific theoretical uncertainties 

The results of Sections 2 and 3 imply that there may be significant theoretical corrections 
beyond NNLO in the standard DGLAP perturbative expansion, which become important at 
either low x, low or low or some combination of these. Indeed there exist several types 
of theoretical correction which may be expected to have such effects. These include 

(i) Higher powers of Inl/a; at higher orders in as, which become important at low x. 

(ii) Increasing powers of ln(l — a;) at higher orders in as, which are generally well understood, 
but which are intrinsically linked to higher-twist corrections. These are important at high 
X and hence low W, see (1). 

(iii) Absorptive corrections which arise from parton recombination and which are higher twist 
in nature. These should become important at low x and Q^. 

(iv) Residual higher-twist contributions, which will be important at low Q^. In practice these 
are often combined with the specific higher- twist contributions, already mentioned, in a 
purely phenomenological parameterization. 

We have performed separate global analyses to investigate the effect of each of these in turn. 
The results are described below. 
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Figure 4: Comparison of the NNLO MRST(cons) partons with the default set MRST(2002) at 
NNLO. The former partons are only reliable for x > 0.005 and > 7 GeV^. 
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4.1 Contribution of higher order \nl/x terms 



It has long been known that the sphtting and coefficient functions typically contain one ad- 
ditional power of \nl/x for each additional power of as- Many of these terms are known ex- 
plicitly [24, 25, 26]. However, a full Inl/x resummation seems to involve many complications; 
for example, treatment of the running coupling, kinematic constraints, coUinear resummations, 
etc. Various procedures exist for incorporating a full Inl/x resummation but there is, as yet, 
no agreed prescription. 

Hence, in the present study, we take a general approach. We include higher order corrections 
to the NLO splitting functions, with the correct maximum power of Inl/x, but we let the 
coefficient be determined by the global fit to the data. We begin by adding one additional term 
to Pgg and to Pqg. For this investigation we choose to include phenomenological ag In^ l/x-type 
terms of the form 



where as = Sas/n and n/ is the number of active quark flavours. Both of the additional terms 
have been constructed so that momentum is still conserved in the evolution. We also add the 
same terms multiplied by Cf/Ca — 4/9 to Pgg and Pqg respectively. This factor of Cf/Ca is 
typical for the results in Inl/x resummation [25]. 

The best NLO global fit, modified as in (2) and (3), gives an improvement in of 21 
compared to MRST2002 [10] and corresponds to A — 3.86 and B — 5.12. These are effective 
coefficients and should not be directly compared with the known coefficients from Inl/x re- 
summations, since they represent the effect of the two towers of In" 1 /x terms. However, the 
values of A and B are of the magnitude expected from the partial information which exists. 
Hence the fit seems to benefit from such terms, which increase the speed of evolution at small 
X. The best-fit value of asiM^) decreases, but only very slightly. The input gluon tiu^ns out 
to be similar to that obtained in the fits when data below x = 0.005 are removed; that is, it is 
larger than the default gluon for large and moderate x, but even more negative at very small x. 
We can only assume that large positive Inl/x terms in the coefficient functions will maintain 
the positivity of observables sensitive to the gluon, such as Fl [27] . Nevertheless, the increased 
gluon evolution at small x does result in a positive gluon more quickly as increases. 

We also investigated the effect of a more flexible parameterization of the \nl/x resummation 
by introducing an additional term, both in (2) and in (3), of the type 






(3) 




(4) 




(5) 
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respectively. Again the form of the terms has been constructed so that momentum is conserved 
in the evolution. With the two additional parameters, we have an effective parameterization of 
both the N^LO and N'^LO corrections at small x. The best fit now has a that is 36 lower 
than MRST2002 with 

^ = -0.27, 5 = 2.79, C = 4.08, D = 1.09. (6) 

The gluon distribution is very similar to the case with just two extra parameters, as is 0:5 (M|). 
However, neither the small x low Q'^ gluon distribution or the parameters A, B, C and D 
are determined very precisely, since there is some trade-off between them. For example, since 
the term in (4) falls off very quickly with (since it behaves as q;|((5^)), a more negative 
input gluon is largely compensated by a larger value of C, resulting in a largely unchanged 
gluon when one moves far away from the input scale Qq — 1 GeV^. However, the generally 
positive input parameters and more negative input gluon is an unambiguous result. The fact 
that ^^(Ml) is not altered much by the inclusion of terms which are important at small x is 
not surprising since it is always the cvohition of the large x structure functions that has the 
dominant influence on the extracted value of 0:5 (M|), and this is not affected by this type 
of modification. The studies with up to four additional parameters were deemed sufficient to 
illustrate the possible effect of In 1/x resummation on the NLO analysis. 

However, we already know that the approximate NNLO splitting functions have a significant 
effect on the small x evolution. These contain one extra power of In l/x in comparison to the 
NLO sphtting functions. Hence we repeated the above analysis at NNLO with the same four 
parameters. A, . . . ,D. The results were largely similar. The improved by 39 from the 
standard MRST(NNLO) fit^ [10], although the parameters take the values 

A = -0.35, S = -2.00, C = 5.49, D = 2.81, (7) 

very different from those at NLO, see (6). In addition they have a markedly reduced positive 
effect for the quarks and a very slightly increased positive effect for the gluon, showing that the 
NNLO contributions (positive for quark evolution at small x but negative for gluon evolution) 
have themselves been important at small x.^ As at NLO, the gluon at high and moderate x is 
increased as compared to the unmodified NNLO fit, while at small x it is more negative. The 
gluons in these ln(l/a;)-modified fits are shown in Fig. 5. At low they are indeed both much 
reduced compared to the default gluons at small x. However, the additional terms in the gluon 
splitting functions and the additional gluons at moderate x drive the small-a; evolution more 
quickly than the default and at high the ln(l/a;)-modified gluons are only a little lower than 
the default gluons at small x. 

Hence, the distinct quality in the improvement in the global fit provides strong evidence 
that large \nl/x contributions beyond even NNLO may be important in the description of the 

^In practice we compared to the NNLO set of partons corresponding to the MRST2002 analysis [10]. 
^The effect is usually presented at quite high Q^, e.g. 30 GeV^, and with partons which are steep at small 

X [20], and is claimed to be quite moderate. At ~ 5 GeV^, with partons which arc flattish at small x, the 
contribution from the NNLO splitting functions is proportionally much larger, and certainly very significant. 
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Figure 5: Comparison of MRST gluons obtained from the fits with additional resummation 
corrections in ln(l/x) with the defauh gluons MRST(2002) at both NLO and NNLO. 



16 



data. However, these corrections are empirically different in the NLO case than in NNLO, 
because the latter has its own large terms for x — > 0. 



4.2 Contribution of higher order ln(l — x) terms 

If we expand the quark coefficient function in powers of as, i.e. 



m=l 



then the coefficient functions c^"*)(a;) and their moments c^'^\N) contain large logarithms as 
X — > 1 {N — > oo) of the form 



ln^-\l -x) 
1 — X 



\k 



' ^^-(In'^A^) (9) 



k 



where k = 1, . . . , 2m. Combining results on soft gluon resummation [28] and finite order results 
[17], Vogt has been able to provide explicit expressions for the first four coefficients [29] in an 
expansion of the form 

oo 
m=l 

(10) 

Hence, we have detailed knowledge of the leading terms in the coefficient functions for the large 
X OT N limit at all orders, and it is argued in [29] that a more efficient convergence is achieved 
if the leading terms are arranged as powers of InA^ rather than ln(l — x). In principle all of 
the terms in (10) should be included. However, the investigation in [29] suggests that, unless 
one is at very high x, or equivalently very low W^, going to a finite order is sufficient. Indeed, 
the suggestion is that for m = 3 the coefficient function is only important above x ~ 0.7, and 
those for m > 3 are only important for x > 0.8. 

In order to investigate this we performed a fit using the NNLO splitting functions and 

coefficient functions, and including also the 0{ag) contribution to the coefficient function in 
(10).^ When using our W"^ cut of 12.5 GeV^, we find that the effect of the 0(a|) coefficient 
function at very high x is almost negligible, since the W"^ cut ensures we are at quite high 
for very large x, and q;|(Q^) is fairly small. In fact the most significant effect of the 0{ag) 
coefficient function is at higher x. Since the coefficient function has a vanishing first moment, 
its positive effect at very high x must be countered by a negative contribution at lower x. In 
practice it increases the structure function for x > 0.55, but decreases it for x < 0.55. This 
decrease is not large, but it affects much more data. In practice the best fit, when including 
the approximate NNNLO coefficient function, is very slightly worse than the usual NNLO fit, 
although the partons and the value of as{Mz) are hardly changed at all. 

^In [29] it is demonstrated that the four terms in (10) are a very good approximation to the full NNNLO 
coefficient function, which can be estimated in detail using the same sort of techniques as in [20]. 
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Prom [29] it is clear that the contributions in (10) beyond 0{ag) are only important at 

even higher x. Hence, we conclude that if we use VT^ut ~ 12.5 GcV^, there is no advantage to 
be gained in including terms beyond the NNLO coefficient function. As we go lower in W"^, 
however, the NNNLO coefficient function does start to have a non-negligible effect. We will 
discuss this in more detail in our analysis of higher-twist corrections. 

Indeed, with reference to higher twist, we should note that there are a number of ambiguities 
when performing large x resummations, not just whether to resum large logarithms in N or 
(1— ,x). We note that the series expansion in (10) is convergent. However, one could alternatively 
calculate the first two towers of terms in the expansion 

dlnaAN,Q') ^ g «™(g2)(^^^ ^ ^ ^^^-i ^ ^ . . ^^^^ 

(^"^^ m=l 

which would give the dominant large- A?" contribution to an effective anomalous dimension for 
structure function evolution. In this case the series would have finite radius of convergence. 
This badly-defined series shows that there are higher twist corrections present, and the ambigu- 
ity in the series is taken as an estimate of the size of the higher-twist corrections in renormalon 
models. The divergence in (11) is reflected in the terms in (10) beyond c„i4 becoming extremely 
large. Strictly speaking one cannot simply perform a large ln(l — x) expansion without encoun- 
tering this interplay with higher-twist corrections, and beyond about NNNLO it is difficult to 
disentangle the two. 



4.3 Absorptive effects 

To investigate the effects of absorption^ we include bilinear terms in evolution equations as 
follows 

dixgix,Q^))_ [.',(.', W (12) 



dlnQ^ lO^R^^'^^^'^'^ ^ 

These terms take into account the dominant small x contributions at lowest order in 0:5, as 
calculated by Mueller and Qiu [30]. We have approximated the two-gluon correlation function 
as 

9^'\x,Q'j) = ^-^^[g{x,Q')r (14) 

as estimated in [30]. We consider two choices of R^, namely R^ = 15 and 5 GeV^^. The 
former represents the naive assumption that R is of the order of the proton radius, whereas 
the latter represents much stronger absorption motivated by 'hotspot' studies which allow a 
large contribution from the diagrams responsible for saturation where both gluon ladders couple 
to the same parton [31]. Of course, inclusion of the absorptive terms leads to a violation of 

^We thank M.G. Ryskin for discussions. 
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momentum conservation. In principle a more complete theoretical treatment would correct for 
this, but in this study we account for the effect by starting the evolution at low scales with 
the partons carrying slightly more than 100% of the proton's momentum. The precise amount 
is chosen so that at high Q^, when the absorptive corrections have completely died away, the 
correct value of 100% is obtained. In practice for i?^ = 15 and 5 GeV~^ we input a total 
momentum of 100.7% and 103% respectively at Qq = 1 GeV^. 

The effect of absorptive corrections in global parton analyses was investigated many years 
ago [32]. The corrections were found to have a sizable, observable effect in the small x region 
accessible to HERA for a parton set B_ in which the gluon was assumed to have a small x 
behaviour of the form xg ~ x~^ at the then input scale of Qo = 4 GeV^. However, with 
the advent of the HERA data, the size of the small x gluon is now known to be considerably 
smaller than that of the B_ set, and hence we may anticipate the shadowing effects will be 
much smaller. 

Here we make two alternative studies. First, we repeat the NLO global analysis, with 
the modifications shown in (12) and (13), starting from MRST2002 partons [10]. For the 

= 15 GeV~^ choice there is almost no change in the fit, the best fit having an essentially 
identical x^- There is no significant change in the partons and as{M^) increases by less than 
0.001. For i?^ = 5 GeV~^ there is a slight deterioration in the best fit of Ax^ — 20. This is 
accompanied by an increase in (M|) to 0.1225 and a slight increase in the small x gluon. 
This shows that the slowing of the evolution by the absorptive corrections is not consistent 
with the data, and the increase in as is necessary to compensate partially for this. 

However, starting from the MRST2002 partons, which have a negative gluon for low x and 
Q^, may be regarded to be inconsistent with an absorptive approach, which really assumes a 
positive input. Hence, as an alternative investigation, we force the input gluon to be positive- 
definite and, indeed, slowly increasing with decreasing x. Since we know that this causes conflict 
with the low data we also raised the cut to Qcut — ^ GeV^. However, after performing 
the fits we investigated the extrapolation to lower Q^. 

In detail, we removed the 'negative' input term [—A_{1 — x)'^- x^^-] from the gluon, and set 
the (conventional) small x power Sg — —0.1, hence ensuring a positive definite starting gluon^. 
For = 15 GeV-2 the best fit, for data above = 5 GeV^, has xg{x, Ql) ~ 0.87x-"-^ at 
small X and is over 200 worse in than MRST2002, most of this deterioration coming from 
the HERA data. The evolution at low x is far too rapid, even though as{M^) reduces to 0.117. 
However, there is also some worsening to the fit to high x data due to the gluon at high x being 
smaller than before. If we extrapolate the fit below = 5 GeV^ then the description of the 
HERA, and even NMC, data becomes very poor indeed, see Fig. 6. 

For = 5 GeV^^ the increased absorptive corrections moderate the problems at lowest x 
and xg{x, Q^) ~ 0.93a;^*^'-'^, i.e., a little larger, and the global is 80 worse than for MRST2002 

^It was also necessary to fix £g to some value, which in practice was chosen to be 0.74, in order to prevent 
the interplay between parameters resulting in an effectively valence-like input gluon distribution. 



19 



MRST(2001) NLO fit , x = 0.00005 - 0.00032 




1 10 

(GeV^) 



Figure 6: The comparison with data of the fit with a fiattish input gluon and saturation 
corrections with the curves extrapolated below the region of validity of the fit. 
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with 0:5 (M|) = 0.118. However, the extrapolation to low is still nearly as bad as in the 
previous case. 

Hence we conclude that if we demand a positive-definite input gluon, which is even slowly 
increasing as x decreases, then absorptive corrections are not sufficient to compensate the faster 
evolution when compared to a negative input (MRST2002) gluon^, and also the necessary 
reduction of the larger x gluon has a detrimental effect on the fit. 

It was hoped that the introduction of absorptive corrections would lessen, or perhaps remove, 
the need for a negative gluon at very low x and low Q^. However we see from the above studies 
that the global fits do not appear to favour the introduction of absorption corrections. This 
result is surprising [34]. We know that at low x about 10% of F2 arises from diffractive events. 
However we cannot simply subtract F2 from the inclusive F2, since the diffractive events which 
originate from low scales are already accounted for in the parameterization of the starting 
distributions at Qq. The diffractive contribution, AF^, arising from higher scales is however 
related to the absorptive correction AF2 to the inclusive F2. The relation is given by the AGK 
cutting rules [35]. When AF2 <S F2, the relation takes the form 



So some absorptive correction is expected to be present. To quantify the amount will require 
an enlarged global analysis incorporating the diffractive structure function data. 

4.4 Higher-twist effects 

In previous papers ([16, 21]) we have examined the effect of including in global fits a simple 
phenomenological parameterization of the higher-twist contribution in the form 



where in practice Di{x) is taken to be a constant, independent of Q^, in each of a number of 
different bins in x. In the fits which include such a parameterization, we have lowered the 
cut to 1.5 GeV^ and the W"^ cut to 4 GeV^. Here we repeat the procedure for our NLO and 
NNLO fits with the most up-to-date data, and also include a fit which has the approximate 
NNNLO coefficient function, as discussed in Section 4.2. This only alters the NNLO results at 
high X. We give results for a LO fit, although we have already seen in [21] that such a fit simply 
fails in many regions of parameter space, and the invoked higher-twist corrections are simply 
mimicking, as best as they can, corrections which really reproduce the NLO (and possibly 
higher) leading-twist contributions. In particular, we see that the large negative higher-twist 
corrections at small x decrease significantly at NLO and even tend to disappear altogether at 
higher orders. The results are summarised in Table 3. 

^This result is in conflict with observations of Ref. [33] in which only LO partons are considered. 



AF2 ~ 




(15) 




(16) 
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X 


LO 


NLO 


NNLO 


NNNLO 


0-0.0005 


-0.07 


-0.02 


-0.02 


-0.03 


0.0005-0.005 


-0.03 


-0.01 


0.03 


0.03 


0.005-0.01 


-0.13 


-0.09 


-0.04 


-0.03 


0.01-0.06 


-0.09 


-0.08 


-0.04 


-0.03 


0.06-0.1 


-0.02 


0.02 


0.03 


0.04 


0.1-0.2 


-0.07 


-0.03 


-0.00 


0.01 


0.2-0.3 


-0.11 


-0.09 


-0.04 


0.00 


0.3-0.4 


-0.06 


-0.13 


-0.06 


-0.01 


0.4-0.5 


0.22 


0.01 


0.07 


0.11 


0.5-0.6 


0.85 


0.40 


0.41 


0.39 


0.6-0.7 


2.6 


1.7 


1.6 


1.4 


0.7-0.8 


7.3 


5.5 


5.1 


4.4 


0.8-0.9 


20.2 


16.7 


16.1 


13.4 



Tabic 3: The values of the higher-twist coefficients Di of (16), in the chosen bins of x, extracted 
from the LO, NLO, NNLO and NNNLO (NNLO with the approximate NNNLO non-singlet 
quark coefRcient function) global fits. 

There are a number of conclusions which can be drawn from the Table. First, we can see 
that there is no clear evidence for any significant higher-twist contributions for x < 0.005. Even 
though it may be argued that the form of the higher-twist corrections at small x is rather more 
complicated than the simple parameterization of (16) (e.g. [36]), they would have to be of 
roughly the same qualitative form, and the lack of any indication of them appears compelling. 

At NLO there is a strong indication of a negative higher twist contribution for x ~ 0.005 — 
0.06. This is required in order to make dF2{x,Q'^)/d\nQ'^ large enough for the NMC data in 
this region. However, the sizes of the D^^s in this region decrease significantly when going to 
NNLO, because both the coefficient functions and sphtting functions at NNLO lead to increased 
evolution in this range, and the evidence for higher twist at NLO seems to be really an indication 
of a lack of important leading-twist, higher-order corrections. 

The main higher-twist corrections appear, as expected, at high x. For x ~ 0.1—0.4 there is 
a shght indication of a negative higher- twist correction at NLO, but this diminishes at NNLO 
and effectively disappears at NNNLO. Hence this is presumably just an indication that leading- 
twist perturbative corrections arc important. A transition is apparent for x ~ 0.4—0.5, and for 
X > 0.5 there is a definite positive higher-twist contribution. However, this contribution has a 
tendency to decrease from one order to the next. Indeed at NNNLO the required higher-twist 
contribution at very high x is similar to that achieved simply from target-mass corrections [37]. 
Moreover, in this very high x domain, the correction to the higher-twist coefficient when going 



22 



from NNLO to NNNLO is as large, if not larger, as when going from NLO to NNLO (the NLO 
to NNLO correction is smaller than the LO to NLO). It is easily verified that this is indeed 
because the NNNLO correction to the structure function is as large at NNNLO as at NNLO for 
this range of x and Q^. This shows that the divergent high-a; perturbative series discussed in 
Section 4.2 reaches its minimum at about NNNLO, and this represents the essential ambiguity 
in this series. Hence the renormalon contribution to the higher twist [38], which comes from this 
ambiguity of the perturbative series, is roughly of the same size as the NNNLO contribution 
(at least in the region of parameter space we are probing). Indeed, the prediction for the 
higher-twist contribution as a function of x in [38] is very similar to that obtained from the 
approximate NNNLO contribution. This implies that at very high x it is pointless to go beyond 
NNLO (or certainly NNNLO) in the perturbative series, since at this order the perturbative 
series and higher-twist corrections become indistinguishable. It is important to realize, however, 
that this is a special feature of high x, i.e., low W"^, and does not imply the same is true for 
other regimes of x and Q^. 

To conclude, by going to higher and higher orders in the leading-twist perturbative expansion 
we see that the only strong evidence for higher-twist corrections is in the region of high x and 
low W^. All shortcomings of low-order perturbative calculations seem to be reduced, if not 
removed, simply by working to higher orders. Moreover, our particular knowledge of high x 
coefficient functions leads us to believe that we have reached the stage where the perturbative 
expansion and higher twist corrections are not really separable. This is further illustrated 
by a method which goes beyond the standard logarithmic resummation using the 'dressed 
gluon exponentiation' approach of [39], which takes into account some all-order information 
on the kernel of Sudakov resummation itself. In Ref. [40], this combined resummation of 
Sudakov logarithms, renormalon contributions and higher twists was confronted by data on the 
Nachtmann moments (i.e. corrected for target mass) of the structure function extracted from 
low data. While it is impossible to estimate precisely the higher twist contributions, one 
allowed description is that where the Sudakov resummation alone is necessary to explain the 
high-moment data. 

5 Uncertainties due to input assumptions 
5.1 Choice of input parameterization 

Over the years we have adopted parameterizations with more and more free parameters, as 
required by the increase in precision and kinematic range of the data, and by the new types of 
data that have become available. Perhaps the biggest single extension of the parameterization 
was the contribution added to the gluon that allowed it to become negative at small x [21]. 
Recently other groups have investigated the influence of parton parameterizations on parton 
uncertainties [6, 8], with various conclusions. In [8], which describes a fit to a fairly limited set 
of data, it is found that there is no real improvement to the quality of the fit once the number 
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of parameters specifying the input parton distributions has increased past a certain number, 10 
in this case. In [6], however, a new type of parton parameterization is used, and it is claimed 
that this is necessary in order to obtain the best fit to the Tevatron jet data. Our studies do 
not support this claim. 

Although the fit to the Tevatron jet data in [6] is indeed better than in [15], and even in 
[10] (where the high x gluon has been improved), we believe very little of this is due to the 
parton parameterization. There are various differences between the CTEQ6 and MRST2002 
approaches to global fits including the choice of cuts, the data sets included in the fit, the 
treatment of errors (particularly those of the E605 Drell-Yan data) and even the definition of 
as{Q^) at NLO. In [41] we discussed in detail the effect of making our starting point for the 
fit more and more like that for CTEQ6, and discovered that when we did so our fit to the jet 
data became of almost similar quality to that of CTEQ6 when using our own parameterization 
(without any of the "kinks" found in some previous best fits to jet data [15]). Some of the 
further cuts we have introduced in the course of the investigations in this paper have led to 
even better fits to the Tevatron jet data. Hence, we do not feel that these particular data 
require us to alter our parton parameters to obtain the best possible fits. However, we do note 
that the fact that our parameterization does allow the input gluon to be negative does have 
important consequences compared to the CTEQ analysis. It means that the CTEQ6 gluon 
is always larger at very small x than ours (even though their starting scale is a little larger, 
Ql — 1.69 GeV^, our gluon is still negative at this Q^), and from the momentum sum rule 
must be smaller elsewhere, in practice at intermediate x. This does lead to predictions which 
are significantly different to ours, examples being the Higgs cross section at the Tevatron and 
LHC as seen in Fig. 15 of [10]. 

Indeed, there is considerable evidence that we have, if anything, too much fiexibility in 
the parton parameterizations. When attempting to diagonalize the error matrix for partons 
when all the parameters were left free we discovered that many were extremely correlated (or 
anti-correlated) leading to some very flat directions in the eigenvalue space [10]. In fact, in 
order to obtain a stable error matrix we found that it was necessary to observe the departures 
away from the best fit while allowing only 15 of our nominal 24 parton parameters to vary. A 
similar effect has, in fact, also been seen in [5], where when producing the error matrix only 
16 out of a possible 22 parameters are allowed to vary, and in [6], where only 20 out of 26 
parameters are allowed to vary. Similar problems are not seen by other groups [7, 8, 9, 3], 
all of whom use smaller data sets and a smaller number of parton parameters, implying that 
they have either the correct number of parameters for the flexibility required by their data, or 
potentially slightly too few.^*^ 

The redundancy observed when calculating the error matrix does not mean that all the 
other parameters can simply be dispensed with, since it is necessary to allow most to vary in 

^''As mentioned above, in [8] the number of required parameters is carefully investigated. In [3], however, one 
can check that more than 2 a variations in the parton parameters would lead to some pathological behaviour 
(particularly for the gluon), implying that the fit is on the verge of the same type of redundancy problems as 
in [10] and [5, 6]. 
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order to obtain the best fit itself, and to allow the partons to have a sufficiently ffexible shape. 
It just means that not all the parameters need to vary in order to investigate small deviations 
from the best fit partons. For example, starting with our previous parameterization for the 
gluon 

xg{x, Ql) = Ag{l - xY^{l + e,a;°-^ + -tgx)x^^ , (17) 
it was necessary to add a term of the form 

-A_{l-xf-x^-, (18) 

in order to let the input gluon be negative at small x. Only 3 parameters were needed to 
investigate small variations, but more parameters are certainly needed in order to obtain the 
gluon for the best fit. However, the total of 7 free parameters we have for the gluon {Ag is fixed 
by the momentum sum rule) do have some genuine degree of redundancy. We have noticed in 
the course of performing many fits that, for any fixed value of eg from —1 to 3, we can obtain 
optimum descriptions of the data which are practically identical in quality. In these fits all 
gluon parameters are significantly different, but the gluon distributions produced are practically 
identical (at least between x — 0.9 and x — 0.00001). This is perhaps the clearest example of 
a single parameter which is essentially redundant, but we have other similar examples. 

Hence, we conclude that our input parameterizations are sufficiently ffexible for the present 
data. We do not seem to have the optimum parameterization for both finding the best fit and 
also investigating fluctuations about this best fit. For a fully global fit, however, no-one else 
seems to have it either. To achieve this would require fewer parameters than at present. This 
might then influence our error analysis using the Hessian approach somewhat, but we feel it is 
unlikely to affect our best fit partons very much at all. 



5.2 Choice of heavy target corrections 

When fitting the CCFR -^2^(2;, Q"^) and F^{x^ Q^) data [42] we have to use some model for 
nuclear shadowing corrections. The form that we use for the heavy target correction factor is 
deduced from a Q^-independent fit to the EMC effect for the scattering of muons on a heavy 
nuclear target (A=56). To be exphcit, we parameterize the correction as 

f 1.238 + 0.203 logio a; for x < 0.0903 
i?HT = < 1.026 for 0.0903 < a; < 0.234 (19) 

[ 0.783 - 0.385 log^o x for 0.234 < x. 

In order to investigate the uncertainty arising from this correction we perform a series of fits, 
maintaining the central plateau in Rwr, but changing the slopes in Inx in the high and low 
X regions. In addition, we allow the normalisation to vary within the experimental error, as 
usual. We study the quality of the fit as a function of the value of the heavy target correction 
at the lowest x value {x = 0.0075), for which CCFR data exist. This data point receives the 
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variation of % with size of heavy target corrections 




0.7 0.75 0.8 0.85 0.9 0.95 



Figure 7: The variation of for global fits with different heavy target corrections to the CCFR 
neutrino data, is plotted against the value of the heavy target correction i?HT at x = 0.0075, 
as explained in the text. The default fit has i?HT = 0.807, whereas for the optimum fit the 
correction is slightly less, i?HT = 0.86. 
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largest heavy target correction, that is i?HT = 0.807 for our standard fit. This is a very similar, 
but larger, correction to that for the highest x value {x — 0.75) for which data exist. Fig. 7 
shows the global ^ function of the value of i?HT at x = 0.0075. Clearly our standard 

shadowing assumption, although near to the minimum, is not the absolutely optimum choice 
of shadowing correction. The data prefers a value of RuTix = 0.0075) of 0.86, i.e. a slightly 
smaller correction than our usual choice, and the improvement in the fit is about 30 units of 
X^. Most of this improvement comes from the fit to the CCFR data on F2^'^\x, Q^), which has 
always had a tendency to lie underneath the theory when shadowing corrections are applied 
(see for example Fig. 9 of [15]). The main effect of the slightly reduced shadowing is just to 
bring these data in line with the theory (although there are some other minor improvements) 
and the parton distributions themselves are not changed much at all, as shown in Fig. 8.^^ 

Similarly the deuterium data is corrected for shadowing effects. These prescriptions are 
not unique. In our fits we normally use the parameterization in [43] which uses a theoretical 
model to estimate the nuclear shadowing corrections at relatively small x. There are other 
alternatives for models of this type of shadowing. There has also been a prescription for 
deuterium shadowing corrections extracted by the SLAC E139/140 experiments [44]. This is 
an empirical extraction which uses a model [45] in which binding effects are assumed to scale 
with nuclear density. This gives a relatively large shadowing correction for deuterium, especially 
at large x, and was used as a basis for an analysis of the d/u ratio in [46]. The validity of this 
extraction is rather controversial (see e.g. [47]), but we use the results simply as an estimate 
on the uncertainty of the deuterium shadowing corrections. As such, a fit with this shadowing 
correction applied to the deuterium data gives an estimate of the model uncertainty from this 
source on the high x valence partons, particularly d{x,Q^)}'^ 

This particular correction leads to a larger neutron structure function at high x. Hence, 
the largest change in the parton distributions is in the high x down quark distribution, which 
increases significantly. As seen in Fig. 9 this increase is a little smaller at high x than the 
uncertainty due to errors on experimental data that was estimated using the Hessian matrix 
method in [10]. From the constraint on the number of valence quarks the increase in dy{x, Q^) 
at high x must be compensated for elsewhere, and indeed we see that it is smaller than the 
default for x = 0.01 — 0.1, but is well within the bounds of the experimental uncertainty. The 
change in the u^{x, Q"^) distribution must lead to a modification of the Uy{x, Q"^) distribution in 
order to obtain the best overall global fit. This modification is shown in Fig. 10. Although it is 
proportionally much smaller than the change in d^{x^ Q^), Uv{x, Q"^) is much better constrained 
by data, and the change due to the different shadowing correction can be as big as, or even 
slightly larger than the uncertainty due to experimental errors on data. However, this is mainly 

"^^The preliminary measurements of F^''^' (x, Q^) by the NuTeV collaboration actually tend to lie a little above 
those of CCFR in the lowest x bins , and hence lie rather closer to the default theory curves with our standard 
shadowing correction [49]. 

^^We note that an examination of theoretical uncertainties due to nuclear effects in the deuteron has recently 
been studied in [48] with the aim of examining isospin depedence of higher twist corrections. 
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Figure 8: The parton distributions obtained from the fit with the optimum shadowing correction 
of Rht{x = 0.0075) of 0.86 compared with the default MRST2002 partons. 
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Figure 9: The ratio of dy{x, Q'^) obtained from the fit with the deuterium shadowing correction 
of [44] compared with the defauh MRST2002 partons, and also with the uncertainty in the 
default d^{x,Q'^) distribution from errors on experimental data found in [10]. 
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Figure 10: The ratio of Uv{x, Q^) obtained from the fit with the deuterium shadowing correction 
of [44] compared with the defauh MRST2002 partons, and also with the uncertainty in the 
default UviXjQ"^) distribution from errors on experimental data found in [10]. 
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so at X < 0.1 where the direct constraint on valence quarks is quite small since sea quarks 
dominate, and the estimated uncertainty from the Hessian method for valence quarks alone 
is likely to have limitations related to parameterizations. Hence, for the down and up quark 
distributions we conclude that model errors on shadowing corrections in deuterium are typically 
of the same size as the errors due to the experimental errors on deuterium structure function 
data. 

Finally we note that the quality of the global fit when using the deuterium corrections 
in [44] is 12 units in x^, better than the default fit. There is a slight improvement in the 
total fit to deuterium data, mainly to the BCDMS data, but a slight deterioration in the fit 
to E605 Drell-Yan data. There is also a slight improvement in the fit to high-£^r Tevatron 
jet data. This is to be expected since the increased down quark distribution at high x leads 
to an increase in the highest Et jet cross-section, which is what is required. However, the 
improvement is limited since an increase in d{x, Q"^) at high x increases the momentum carried 
by the down quark. From the momentum sum rule this makes it harder to have a large gluon 
distribution at high x, which the fit would like. In practice a compromise is reached, i.e. an 
even greater enhancement in d^{x^Q'^) at high x would slightly improve the fit to deuterium, 
and so other data, but actually makes the fit to Tevatron jet data worse. The relatively small 
improvement in the global from this model of deuterium shadowing, partially limited by 
tension between different data sets and parton distributions, leads us to conclude that there is 
no strong supporting evidence for the model. Nevertheless, the small pull of the evidence is in 
favour of some deuterium shadowing at high x. 

5.3 Size of input strange sea 

Global fits have traditionally assumed that the shape of the input strange quark sea distribution 
is the same as the average of the input u + d distribution. The primary experimental constraint 
on the strange distribution is provided by data on dimuon production in neutrino-nuclei deep 
inelastic scattering [50, 51]. These data are consistent with the shape assumption, and in 
addition constrain the magnitude of the strange distribution to be approximately half of {u + 
J)/2 at low Q"^. For this reason, we conventionally choose 



and performing global fits, including the 'data' on the strange sea provided by CCFR [50], for 
various different fixed values of the parameter k. The variation in is shown as a function 
of K in Fig. 11. The continuous and dashed curves correspond to including all CCFR data or 



s{x) = -{u{x) + d{x)) 

at Qo = 1 GeV^. We now investigate the sensitivity to this input assumption. 
First, we check the validity of our choice by expressing 



(20) 




(21) 
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Figure 11: The variation of for global fits witfi different values of k, defined in (21), which 
gives the strength of the strange quark sea relative to the non-strange sea. The two curves with 
the deep minima include the CCFR dimuon data [50] in the global x^, whereas when these 
data are omitted the shallower x^ profile is obtained. 
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omitting those CCFR data at the lowest value, Q"^ — 1 GeV^. The latter may be more 
appropriate since NLO leading-twist distributions might become unreliable at such a low scale. 
We see that our default choice oi k = 0.5 is near the minimum when the = 1 GeV^ CCFR 
data are included, but is perhaps a little large when these data are omitted. The choice k = 0.44 
appears optimum for both data selections. We therefore make a set of partons available with 
this value. In the future, NuTeV dimuon data will better determine the strange sea. 

We also investigate the effect of the variation of n on the quality of our standard global fit 
(in which the contribution of the dimuon data is not included). This is shown by the (shallower) 
dash-dotted curve in Fig. 11. We see that the minimum is obtained for k = 0.38 and that the 

is about 8 lower than for the MRST(2002) set. Hence this is a rather small improvement 
in the fit and leads to a very small change in the parton distributions. This value of k is near 
the limit of acceptability as determined by the dimuon data. We conclude that the level of 
uncertainty of parton distributions, and related quantities, is rather small, particularly as it 
weakened by evolution to higher Q^. 



5.4 Possible isospin violations and s s 

Isospin symmetry implies that the parton distributions of a neutron are obtained from those of 
the proton simply by swapping the up and down quark distributions, i.e. d"'{x, Q"^) = u^{x, Q^) 
and m"(x, Q^) = dP{x,Q'^). In the absence of any obvious evidence to the contrary this is 
always assumed to be true in global fits. There are many sets of data in the global fit which 
would in principle be sensitive to any isospin violation. These are the various sets of deuterium 
structure function data (SLAC, BCDMS, NMC), the CCFR neutrino structure function data 
from isoscalar targets, the E605 Drell-Yan data on a copper target, and the NA51 and E866 
Drell-Yan asymmetry^^ data. However, in the global fit, all of these data sets are well described, 
implying that isospin symmetry breaking is small. This may be illustrated as follows. Assuming 
that the structure functions are dominated by the up and down quarks (and antiquarks), which 
is largely true, then the measurements of the 7-exchange contributions to F2 determine the 
quark contributions 

oc - (m" + m") + + r) (23) 
9 9 ^ ' 

whilst the structure functions for u and p interactions on an isoscalar target constrain the 

combinations 

o^dF + vF + <r + vJ' (24) 

F^ oc + + m" + J" (25) 

xF^ (X(F -vP ^cT -vJ" (26) 

^^The pp and pn asymmetry measurements involve Drell-Yan production on deuterium, as well as hydrogen, 
targets. 
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xF^ (X uP - dP + u"" - (T. (27) 

However the CCFR data [42] for the neutrino structure functions F2 and xF^ that are used in 
the global analyses average over the v and P interactions^^. Note that the F^jF^ NMC data 
determine the ratio of the quark combinations of (23) and (22) with about a 2% error. The 
data thus allow little flexibility in the neutron parton distributions even without the additional 
constraints of Drell-Yan data and VF-asymmetry data (which constrain the proton valence and 
sea quarks). 

We consider the possibility of isospin violation in the valence quarks and the sea quarks 
separately. For the sea quarks we assume 

<ea(^)=Ca(^)(l + '^), (28) 
Ca(^)=<a(^)(l-<^)- (29) 

This type of violation is expected from theoretical models, and is consistent with momentum 
conservation, up to very small violations due to a non-zero value of {d^f.^ — 'uP^^^- Strictly 
speaking it is not preserved by evolution of the partons, but in the kinematic regions of interest 
the violation is very small. Somewhat surprisingly, we find that a certain amount of isospin 
violation of this type is actually preferred by the data. The best fit is obtained for b = 0.08, 
i.e. an 8% violation of isospin in the sea. This fit has a total 20 better than the default 
best fit at NLO. The majority of this improvement comes from the fit to NMC data on F^f /Ff , 
which is raised a little by the increase in Ug^^{x). The fit to the E605 Drell-Yan data is also 
markedly improved. The change in the up and down sea quarks for this fit compared to our 
default partons is shown in Fig. 12. This clearly shows the preference of the data for the up 
sea distribution in the neutron to be enhanced, and the down quark suppressed. 

An increase in of 50, corresponding to a 90% confidence limit from our arguments in 
[10], arises when 5 — 0.18 or 6 — —0.08, see Fig. 13. In the former case the NMC data on 
F2/F2 (which is the most sensitive discriminator of isospin violation) has now become too 
large, the fit to BCDMS data has deteriorated^^ and the description of the E866 Drell-Yan 
asymmetry has become very poor. The fit to E605 data has continued to improve, however. 
For 6 = —0.08 the prediction for F2/F2 has become too small, and the fit to E605 Drell-Yan 
data has deteriorated. 

For the valence quarks a similar model of isospin violation is not possible because it would 
violate the valence quark number counting for the neutron. Hence we consider a violation of 
the type 

<ix) = dl{x) + Kf{x), (30) 

^^When available, the NuTeV measurements of the v and 9 interactions individually will offer even more 
stringent checks of the isospin relations. 

^^The deterioration is caused by the decrease of dl^^{x), which arises in order to prevent u^^^{x) becoming 
too large. 
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Figure 12: The up and down sea quarks for the best fit with allowed isospin violation in the sea 
quarks compared to the default partons. The ratio is shown both for the proton and neutron 
sea quark distributions. 
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where f{x) is a function which has zero first moment. A suitable function, which has the same 
type of behaviour as the valence quarks at high and low x is f{x) = (1 — x)^x~^'^(x — 0.0909). 
However, there is a further constraint. If we add Kf{x) in order to obtain u^{x), then we must 
subtract nf{x) in order to obtain d^{x), i.e. 

<(x)=<(x)-«/(x), (31) 

in order to ensure momentum conservation. Hence, we consider isospin violation of this form. 
Again, the isospin violation is not exactly preserved by evolution, but is correct to a very good 
approximation . 

For valence quarks, there is very little preference for isospin violation. The best fit is 
obtained for k — —0.2, see Fig. 13, but this gives an improvement in of only 4. It corresponds 
to at maximum about a 3% violation of isospin for u^{x). The 90% confidence level is obtained 
for K — —0.8 or K — 0.65. In the former case, the main deterioration is in the description of the 
CCFR F2{x, Q^) data in the region of x = 0.2, and in the latter case, it is CCFR F!^{x, Q^) data 
and BCDMS F^{x, Q^) data which are both badly fit at X ~ 0.5. For positive k, d^{x) decreases 
at high X to compensate the isospin violating term and in order to fit the ^^(a;, Q^) data u^^x) 
decreases, but less severely (due to the higher charge weighting). Hence, the failure in the 
fit occurs when u'^{x) has increased too much for the BCDMS deuterium data but the larger 
decrease in d^{x), compared to the increase in u^{x), leaves the prediction for F^{x, Q^) (which 
is sensitive to u^{x) -\-d^{x)) too small. For negative d^{x) increases at high a;, to compensate 
the isospin violating term, and in order to fit the F^{x, Q'^) data, uP{x) increases less severely. 
In this case u^{x) + d'^{x) decreases, and is too small for F2^(x, Q^). These upper and lower 
limits represent an isospin violation of at most C(10%) for m"(x). We offer no theoretical model 
for why our isospin violation is of the form seen, and note that calculations of the effect, e.g. 
[52], often indicate smaller results. However, these calculations are very difficult, depending 
on intrinsically nonperturbative physics, and relying on models and assumptions. We simply 
examine the empirical evidence given by the data. 

It is interesting to compare the level of isospin violation with that needed to explain the 
NuTeV sin^ 9w anomaly [53] . The quantity measured^^ by NuTeV is 

R- = ^NCJl^c. (32) 

^CC ^CC 

In the simplest approximation, i.e. assuming an isoscalar target, no isospin violation and equal 
strange and anti-strange distributions, this ratio is given by 

R" -sin^Ow, (33) 

and so the measurement gives a determination of sin^^i^^. NuTeV find sin^^^ = 0.2277 ± 
0.0013(stat.) ± 0.0009(syst.) [53], compared to the global average of 0.2227 ± 0.0004, i.e. about 

^^The NuTeV experiment does not exactly measure in part because it is not possible experimentally to 
measure neutral current reactions down to zero recoil energy, see [53, 54] . 
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a 3(7 discrepancy. However, if one allows for isospin violation then the simple expression becomes 
modified to 

sin^ % + (1 - ^ sin^ (34) 

where 

[SU^] = [' x«(x) - <(x)), [SD^] = [' x«(x) - <(x)), (35) 

Jo Jo 

are measures of the inequality in momentum fraction of the valence quarks induced by isospin 

violation, and [V~] ~ 0.45 is the overall momentum fraction carried by the valence quarks. One 

can easily see that given a fixed value of measured R~ , a negative value of k, moves the extracted 

value of sin^ 9w downwards. The approximate effect can be found simply by using Eq. (34), but 

a more precise result is obtained by applying the functionals presented in [54] , which account for 

the complications of the measurement. This reduces the naive modification by ~ 10%. Using 

our best fit value oi k, — —0.2 we obtain [SUy] — —[SDy] — 0.002 (corresponding to < 0.5% of 

the momentum carried by the valence quarks) and the modification is Asin^Ow = —0.0018. 

Hence, about la — 1.5a of the 3a discrepancy is removed. The determination of sin^ 9w is far 

less sensitive to the isospin violation in the sea quarks. Our preferred violation with 6 = 0.08 

is in the wrong direction to account for the discrepancy, but only reduces the effect from the 

valence quarks slightly to Asin^ 9w = —0.0015. Hence, the total result from the best fits with 

allowed isospin violation is to reduce the 3a discrepancy to a 2a discrepancy. However, the 

allowed range of isospin violation could easily allow the discrepancy to be removed altogether, 

or even to be made worse. It is nevertheless interesting that the weak indication given by the 

data in the global fit is such as to reduce the discrepancy by a significant amount. 

One can also investigate the possibility of s{x) ^ s{x). The only data in the fit which are 
sensitive to this difference are the NuTeV dimuon data [51]. These data are more difficult to 
analyse than the similar CCFR data because they are presented in a much more exclusive form. 
However, in principle such data will allow for a much more detailed analysis. The NuTeV group 
themselves have performed an analysis of their data allowing the six) and s(x) distributions 
to have different normalisations and different (1 — x)^ behaviour at high x. Their analysis 
indicates that the data would prefer a slight (11%) excess of s{x) over s{x) [51]. However the 
analysis is at leading order, and does not impose the quark counting rule, i.e. equal number of 
strange and anti-strange quarks. It is complicated to improve this analysis to the full NLO level. 
CTEQ have performed similar preliminary analyses [55], obtaining very similar results. Most 
predictions of physical quantities are insensitive to the potential imbalance of s{x) and s{x\ 
but the NuTeV sin^ anomaly is affected by any difference [53]. However, the excess of s{x) 
over six) actually makes the discrepancy slightly worse. More precise data and a full theoretical 
treatment will hopefully lead to an improved understanding of this question in future. 

^^CTEQ have produced fits where the quark counting is imposed and find that the excess in over s(x) 
in the region of data must then be countered by an excess of six) over six) at higher x, leading to a positive 
momentum excess of s(a;) over six) [56]. This is then in the correct direction to reduce the NuTeV sin^^^^ 
anomaly. 
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Figure 13: profiles showing the effect of isospin violation in the sea and valence quark sectors 
respectively. 6 and k are defined in eqs. (28) and (30) respectively. 
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Ois(Mz) ± expt ± theory ± model 

NLO" 

CTEQ6 
ZEUS 
MRST03 
HI 

Alekhin 

MRST03 5 0.1153 ±0.002 ±0.003 
Alekhin 1 0.1143 ± 0.0014 ± 0.0009 

Table 4: The values of as{M^) found in NLO and NNLO fits to DIS data. The experimental 

errors quoted correspond to an increase Ax^ from the best fit value of x^- CTEQ6 [6] and 
MRST03 are global fits, where the latter correspond to the 'conservative' sets of partons of 
Sections 2.5 and 3.2. HI [7] fit only a subset of data, while Alekhin [57] also includes Fl*^ 
and ZEUS [9] in addition include xF^ data. 

6 Implications for as and predictions for observables 
6.1 Determination of as 

As we have already demonstrated, there is a significant amount of variation in our extracted 
value of ^^(Ml) when we vary the input assumptions for our fitting procedure, particularly 
when we vary the x and cuts. Previously we have always determined «5(M|) from a global 
fit using our default cuts, that is no cut on x and a cut of 2 GcV^. As demonstrated earlier, 
the evidence suggests that when fitting this full range of data, standard NLO (or NNLO) 
perturbation theory is not completely sufficient, and the variation in as{M'^) is due to the 
parameters in the fit compensating for the deficiencies of the theoretical treatment. Hence, the 
'true' value of as{M^) should be that corresponding to the 'conservative' partons at both NLO 
and NNLO. We present these values below in Table 4, labelled MRST03, together with other 
recent determinations. Note that the conservative NLO partons themselves use 0:5 = 0.1162. 
However the profile versus 0:5 is very fiat between 0.116 and 0.117. We therefore choose the 
mid-point as the best value of as, and Ax^ = 5 to give the la error of ±0.002. 

As mentioned in Section 2, a Bayesian approach to determining parton uncertainties [2] 
gives as{M'^) = 0.112 ± 0.001. However, in order to satisfy the strict requirements of consis- 
tency between the data sets, this Bayesian analysis only uses the BCDMS [13], E665 [14] and 
Hl(94) [12] data for Ff. Bearing in mind that the Hl(94) data have relatively large errors. 



100 0.1165 ±0.0065 

50 0.1166 ±0.0049 ±0.0018 

5 0.1165 ±0.002 ±0.003 

1 0.115 ± 0.0017 ±0.005lgS 

1 0.1171 ±0.0015 ±0.0033 
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the value of as{M]r) that is obtained simply reflects the original BCDMS determination of 
Ref. [62]. 

Prom Table 4 we see that the various determinations of 0:5 (M|) have approximately con- 
verged to a common value, even though they are based on different selections of the DIS and 
related data. Averaging the two global NLO analyses we have 

asiMl) = 0.1165 =b 0.004. (36) 

Previously, the MRST value [15] was larger, as{Mz) — 0.119. This was due to the attempt to 
fit data in regions where the theoretical corrections to NLO DGLAP appear to be important. 
It is well illustrated by Fig. 16 of [15], where we see that the optimum value of ^^(Ml) varies 
considerably from data set to data set. However, when the 'conservative' data cuts are applied, 
the tension between the data sets is reduced enormously, as can be seen by comparing Fig. 14 
with Fig. 16 of [15]. Some data sets, particularly the SLAG data, which prefer a high value 
of as{Mz), and the E605 Drell-Yan data and the BCDMS data, which prefer a low value of 
as{Mz), still pull strongly away from the minimum value. However, the other data sets are 
now at, or near, their minimum for the best fit value of as{M^), which was certainly not 
the case for the default fit. The DO and CDF jet data are a particularly good example, where 
not only are the data better fit by the 'conservative' sets than the default sets, but they are 
no longer pulled to extreme values of as{M^) in order to obtain their best individual fits. The 
same marked improvement occurs for the NMC data which survive the x and cuts. This 
increased compatibility between the data sets, and also between the data and the theory, is why 
the tolerance Ax^ has been reduced from 20 (in [15]) to 5 in the present study. The reduction 
in tension between data sets also occurs when ln(l/a;), ln(l — a;) or higher- twist corrections are 
included, as discussed in Section 4. 

There are fewer extractions of the value of as{Mz) using NNLO global or semi-global fits to 
DIS and related data. The results are also shown in Table 4. Again we see good agreement, and 
a small, but definite, reduction from the NLO value. We found, in this case, far less sensitivity 
to the data cuts, indicating that some important theoretical corrections are already accounted 
for at NNLO. 

6.2 Predictions for W and Higgs hadroproduction 

Predictions for physical quantities are, like the value of 0:5 (M|), sensitive to the 'theoretical' 
uncertainties in the global parton analysis. For illustration, we show in Figs. 15 - 18 the 
variation of the W and Higgs cross section predictions for the Tevatron and LHC as a function 

of Xcnt and Ql^^. 

The cross sections for W and H production at the Tevatron sample partons down to x ~ 
0.005, and so are only directly sensitive to partons within the range of our most conservative 
cuts. However the cross sections can still vary if we change the values of the cuts due to the 
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Figure 14: The quality of the fit to the individual data sets included in the NLO global analysis 
with Xcut = 0.005 and Qcut — 10 GeV^, shown together with the overall total function 

of «5(M|). 
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Figure 15: Predictions for the {Mh = 120 GeV) Higgs cross section at the Tevatron {a/s — 
1.96 TeV) at NLO and NNLO for various values of Xcut, and for the 'conservative' partons with 
a cut on both x and (shown as open symbols). 



pq 



2.80 
2.75 
2.70 
2.65 



2.60 - 



2.55 - 



2.50 



MRST NLO and NNLO partons ^ @ TcvatrOn - 



NNLO 



NLO 



♦ ♦ 

O Q'cut = 7GeV' 

• • 

O Ql =10GeV^ 




Figure 16: Predictions for the W cross section (times the leptonic branching ratio Bi^ = 0.1068) 
at the Tevatron {^/s = 1.96 TeV) at NLO and NNLO for various values of Xcut, and for the 
'conservative' partons with a cut on both x and (shown as open symbols). 
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Figure 17: The same as Fig. 15, but for the LHC energy of = 14 TeV. 
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readjustment of partons above Xcut a-nd <5cut' Higgs production, due to the different 

values of 0:5 (M|) extracted. This is evident in Figs. 15 and 16. 

The NLO prediction^^ for the cross section for the production of a Higgs of mass 120 GeV 
at the Tevatron rises steadily as Xcnt is increased, saturating at an increase of about 9% for 
Xcut — 0.005. This rise is due to the increase of the gluon distribution at moderate x values as 
the value of Xcut is raised, see Fig. 1. There is a slight decrease in the value of ^^(Ml) with 
increasing Xcut- but the effect of this on the Higgs cross section is completely outweighed by 
the large increase in the gluon in the relevant x range. We also sec from Fig. 17 that, when 
Qcut is increased to 10 GeV^, corresponding to our conservative set of NLO partons, the Higgs 
cross section is only increased by 1%. This is mainly due to the drop in as{Mz) to 0.116, but 
also due to a slight decrease in the gluon for the relevant value of x ~ 0.06 (for central rapidity 
production), compared to the case where only the cut in x is applied. This arises because the 
g2 cut of 10 GeV2 ehminates the NMC F^^'^\x,Q'^) data for X ~ 0.05, which prefer a steeper 
rise of F2{x,Q'^) with Q^, and hence a larger gluon in this range (as well as a larger value of 
as{M^)). Nevertheless the value for the 'conservative' set represents a non-negligible increase 
in the Higgs cross section compared to the prediction of the default set. 

From Fig. 16 we also see that the NLO prediction for the cross section for W production at 
the Tevatron also rises steadily as Xcut is increased, saturating at an increase of about 2% for 
Xcut — 0.005. This is due to the increased evolution of the quarks driven by the increase of the 
gluon in the relevant range. Again, when we also impose Qcut — i-0 GeV^ the predicted value 
of the cross section decreases. The increased value of Qcut allows the input quarks to be larger 
for X ~ 0.05 (too large for the NMC data now cut out ), and the improvement in the quality 
of the fit requires less increase in dF2/dlnQ'^. Compared to the quarks at x ~ 0.05 obtained 
using Qcut = 2 GeV^, the quarks corresponding to Qcut = 10 GcV^ cross in magnitude in the 
region of = 250 GeV^ and become 0.4% smaller at ^ m^. 

At NNLO we see from Figs. 15 and 16 that the predictions for W and H production at the 
Tevatron are much more stable to variations of x^ut- The reason is that the increased evolution 
of quarks for x ~ 0.05, due to the NNLO contributions to the splitting and coefficient functions, 
requires much less change in the gluon, and consequently in dF2/dlnQ'^ in this range. The 
additional imposition of Qcut — GeV^, however, reduces the Higgs cross section significantly, 
by about 5%, even though at NNLO it leads to only a slight decrease in the value of as'(M|). 
This is because the loss of the NMC F^^'^^x^Q^) data for x ~ 0.05 below Q"^ of 10 GeV^ 
has resulted in a large reduction in the gluon for x ~ 0.1, as seen in Fig. 4. The prediction 
for the W cross section also falls, by about 1%, for the same reason as at NLO. Hence at 
NNLO the 'conservative' parton set predicts a small decrease in both an and aw compared 
to the prediction of the default parton set, while at NLO there is a small increase in both. 

^*For the Higgs cross section calculations described in this section we use the full NLO QCD correction in the 
nit ^ Mh limit [58] and the soft-virtual-coUinear x-spacc (SVC a;) approximation to the full NNLO correction, 
again for toj 3> Mh, taken from Ref. [59]. In both cases the factorisation and renormalisation scales are set 
equal to Mh- 
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As a consequence, the 'conservative' parton predictions show shghtly greater convergence with 
increased perturbative order, than the default predictions. 

At the LHC, the cross sections for W and H production at central rapidity sample partons 
a,t X — 0.006 and x — 0.0085 respectively, so these are safely predicted using the 'conservative' 
partons. However, the total cross sections sample a wide range of rapidity, and in fact probe 
down to a; ~ 0.00008 for W production and x 0.0006 for H production, and hence the 
predictions using the 'conservative' partons arc not guaranteed to be reliable. Wc present such 
predictions to exhibit where the reliability docs break down, and to what extent. 

We notice (Fig. 17) that the predictions for the Higgs cross section at the LHC are actually 
rather stable to changes in the value of x^ut at both NLO and NNLO. From Fig. 3 this is not 
too surprising. At the high scales that are needed for Higgs production, the NNLO gluon is 
hardly changed compared to the default. The NLO gluon, with Xcut = 0.005 applied, is a little 
larger in the central rapidity region, and only falls to a much smaller value than the default 
at values of x that are only making very small contributions to the total cross section. Hence, 
there is only a shght increase in the prediction for the Higgs cross section compared to the 
default. When the additional cut in is applied, in order to obtain the conservative sets, 
both the NLO and NNLO predictions decrease slightly due to the decreases in as^M^)- In 
fact, at NLO and NNLO the predictions using the 'conservative' sets finish very close to the 
default predictions (-1.5% and -1% for the changes at NLO and NNLO respectively). This 
implies that there is little theoretical error associated with Higgs production at the LHC due 
to the partons themselves. Indeed the variation with x^ut and Qcut LHC, exhibited in 

the figure, is evidently much smaller than the correction in going from NLO to NNLO, that is, 
smaller than from the NNLO coefficient function contribution.^^ 

For W production at the LHC, Fig. 18, the story is rather different. At NLO there is a 
steady, and rather dramatic decrease in the predicted cross section as Xcvx is increased. This 
culminates in a drop of 20% for Xc^t — 0.005 (which is insensitive to the additional cut). 
The reason for this is clear from examination of Fig. 2. For the sets with Xcut = 0.005 the quark 
distributions for x < 0.005 are much smaller than those of the default set. Also the gluon at 
small X and low is reduced, and therefore the evolution of the quarks is slower than the 
default set. Hence, although at ~ 10*^ GeV^ the quarks at the central rapidity value of 
X = 0.006 are actually much the same as in the default set, they quickly become reduced at 
smaller x, and very much so for x < 0.001. However, much of the total cross section comes 
from W rapidities corresponding to such low x for one of the two quarks contributing to the 
reaction, and, as a consequence, the contribution to the cross section at high rapidities is much 
reduced. This is illustrated in Fig. 19, which shows the differential cross section as a function 
of the W rapidity. At central rapidity, the cross section is even increased slightly, but it falls 
away very quickly with increasing \yw\, resulting in the 20% loss for the total. However, at 
NNLO there is far greater stability. There is a slight decrease in predicted cross section with 

^^By way of calibration, the scale dependence of the Mh = 120 GcV Higgs cross section at the LHC is quoted 
as ±10% (at NNLO) and ±8% (using NNLL soft-gluon resummation) in the recent study of Ref. [60]. 
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Figure 19: The predictions for the rapidity distribution of the W cross section at the LHC for 
both the default MRST2002 set and for the 'conservative' set. 
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increasing Xcut, which can be understood from the small decrease in the quarks for x in the 
region 0.0001-0.005, as compared to the default set, as seen in Fig. 4. However, as we have 
already commented for the gluon, the change in the very low x partons with increasing Xcut is 
very much reduced at NNLO, because the NNLO contributions themselves include important 
corrections to the small- a; behaviour of the partons and structure functions. 

We believe that the variation of the predicted cross sections with the value of Xcnt and Q^ut 
gives a rough indication of the theoretical uncertainty due to the parton distributions. From 
the NNLO prediction for the cross section for W production at the Tcvatron, we sec from 
Fig. 16 that, not only is the value stable, but that the theoretical uncertainty is of order 1-2%. 
This is similar in magnitude to the uncertainty due to experimental errors, which was obtained 
in Ref. [10]. This gives a total error of about ±2%, which would imply that observing the 
W production rate at the Tevatron (and comparing with the NNLO prediction) could serve 
as a valuable luminosity monitor. The theoretical error, from the partons, for NNLO Higgs 
production at the Tevatron is rather larger (about 5%), due to the adjustment of the gluon 
over the whole x range and changes in as, when the input assumptions to the fit are varied. 

The variation at the LHC is relatively small for the Higgs cross section, which does not 
probe partons at too low x, and is swamped by the higher-order corrections to the partonic 
cross section. However, it is potentially much larger for the W cross section at the LHC, which 
relies on lower x partons. The 20% variation in the W cross section at the LHC at NLO 
implies that important theoretical corrections are required. The reduction of this variation 
to ~ 3% at NNLO implies than much of this correction has occurred in going from NLO to 
NNLO. However, in this case, the uncertainty is still larger than the change due to the NNLO 
contribution to the partonic cross section and is the dominant theoretical uncertainty. ^° The 
results in Section 4 imply that even at NNLO additional theoretical precision could be obtained 
by further theoretical corrections, e.g. a correct inclusion of ln(l/a;) terms at higher orders. 
However, NNLO seems to be a great improvement on NLO, if one wishes to predict quantities 
sensitive to partons for x much lower than 0.005. 



7 Conclusions 

In a previous paper [10] we have already studied the uncertainties of parton distributions, 
and related observables, arising from the errors on the experimental data used in the global 
parton analysis. However in that paper we had already commented that in many cases the 
major uncertainty could be due to corrections to the standard DGLAP evolution and due to 

^''A very recent prediction of the W and Z cross sections at NLO and NNLO, with errors, has appeared in 
[61]. This disagrees with our predictions by 0(5%), which is larger than the total quoted error in [61] (2% for 
the Tevatron and 3% for the LHC). The main reason for this discrepancy is almost certainly due to the absence 
of a number of sets of data which determine the precise form of the quark distributions in the partonic fit [57] 
that is used. Hence, in our spirit of determining parton distributions, we would deem this to be a removable 
discrepancy. 
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assumptions used in the fit procedure, in other words due, collectively, to so-called theoretical 
errors. In this paper we have studied a wide range of sources of theoretical error and their 
potential consequences. 

To begin, we investigated the possible corrections to standard fixed-order DGLAP analyses. 
The investigation was performed in two alternative ways. First we made an empirical study in 
order to find those kinematic regions where DGLAP evolution was fully consistent. This was 
done, at both NLO and NNLO, by gradually eliminating data until the analyses were stable 
to further cuts. We found that this was best achieved by imposing separate cuts on x, Q'^ and 
of the data fitted. 

• For W'^, stabihty was achieved by raising the cut from our default value of 12.5 GeV^ 
to just 15 GeV^. However, in the region of low W'^, we noted that there exists some 
incompatibility between the data sets. 

• For X, stability was achieved only for the relatively high value of Xcnt = 0.005. When 
this cut was applied, the gluon distribution above this x value increased, relative to the 
default set, at the expense of the loss of the gluon at smaller x. This improved the quality 
of the fit to the Tevatron jet data, and also to the structure function data for x ~ 0.01 
due to an increase in dF2/d\nQ'^. The Xcut required was the same at NNLO as at NLO, 
but the modification of the gluon was considerably smaller at NNLO. 

• For Q^, stabihty is reached for Q^ut ~ 7-10 GeV^ at both NLO and NNLO, the conver- 
gence to stability being quite gradual as Q^ut raised. The slow convergence indicates 
that higher-order corrections, rather than higher-twist effects, are important at relatively 
low Q2. 

At both NLO and NNLO we also considered combinations of the above types of cuts. 
Thus we found the full kinematic region where fixed-order DGLAP analysis is appropriate, 
together with the corresponding sets of conservative partons. At NLO the domain is given by 

> 15 GeV^, Q2 > 10 GeV^ and x > 0.005, whereas at NNLO it is given by > 15 GeV^, 

> 7 GeV^ and x > 0.005.^^ Within these regions these conservative sets of partons are most 
reliable, but should not be used outside the domain. Indeed, outside the region they may be 
completely incompatible with the data. However, we note that, although the NLO and NNLO 
regions of stabihty are similar, the NNLO partons are far closer to the default partons outside 
the stable domain, indicating smaller theoretical uncertainty at NNLO. 

Gomplementary to the above empirical study, we also made an explicit investigation of the 
following variety of possible theoretical corrections. 

^^In practice the conservative partons were obtained using W|?„^ = 12.5 GeV^, but raising Wj?^^ to 15 GeV^ 
has a neghgible effect on the partons. 
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• Higher-order parametric ln(l/a:) corrections to the sphtting functions were found to im- 
prove the quahty of the fit. This occurs both at small as expected, and also at larger 
X where the partons are allowed to re-adjust to a form similar to that of the conser- 
vative parton set (with Xcut > 0.005). The required corrections were not the same at 
NNLO as at NLO, being overall a little smaller at NNLO, reflecting the important small 
X contributions introduced at NNLO. 

• Shadowing (or absorptive) corrections were introduced, but found to have little effect 
when the default partons were used as the starting point. If the input gluon was forced 
to be even very slowly increasing at small x, then evolution is too rapid even with shad- 
owing corrections included. Further study would require the simultaneous description of 
diffractive structure function data. 

• Parametric higher-twist contributions were introduced and quantified. These decreased 
as we progressed to higher perturbative order, and, at NNLO, were only evident for 
X > 0.5. However, at both NLO and NNLO, the change in the input partons is similar 
to that invoked by increasing Qcut- The high ,t, or equivalently low W"^, higher-twist 
contribution is intrinsically related to the ln(l — x) resummation. This is because the latter 
is inherently divergent and the ambiguity in the perturbation series must cancel that in 
the higher-twist contribution. We found that the effect of the well-determined NNNLO 
high X contribution was important, and reduced the high x, higher-twist contribution 
still further. However this NNNLO contribution demarks the perturbative order beyond 
which the series fails to converge. Hence we conclude that resummation beyond NNNLO 
becomes indistinguishable from higher-twist corrections at high x. 

As well as the theoretical corrections to the standard DGLAP evolution there arc other 
potential sources of uncertainty in the fitting procedure. Among these we considered the 
following. 

• We found our input parameterization was sufficiently flexible to accommodate the data, 
and indeed there is a certain redundancy evident. Hence we conclude that the form of our 
parameterization does not significantly constrain the description, though a more efficient 
parameterization may be possible. Allowing our input gluon to be negative at small x 
does have important consequences however. 

• Heavy-target corrections are required when fitting to neutrino data. We found that, 
although our default choice was not quite optimum, changes resulted in different data 
sets becoming more compatible with each other, and led to only minimal changes in the 
partons. 

• Shadowing corrections are also required when fitting to deuterium structure functions. 
We found that there was a small amount of evidence for some high x shadowing, and that 
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the uncertainty in valence partons, particularly dy{x, Q^) at high x, due to the model 
uncertainty in deuterium shadowing is similar to the uncertainty due to the experimental 
errors on the data. 

• The input strange quark sea distribution is parameterized as K{u + d) /2, where the default 
choice was k, — 0.5. We found that both the global fit and the description of the CCFR 
dimuon data prefer a smaller value of k ~ 0.44. This choice, which gives a smaller strange 
sea distribution and slightly larger u, d, will be implemented in future fits. However it 
leads to only very small changes in most physical quantities. 

• The possibility of isospin violations was investigated in both the sea and valence quark 
sectors. In the former case, we found an that increase of u'^^^^ (and a corresponding decrease 
of (i"ea)) comparison to the default set, was preferred by the data at the 8% level. In 
fact the acceptable range of u^^^ allowed an increase up to 18% and a decrease down to 
8%. For the valence quarks there was no significant improvement in the fit due to possible 
isospin violations, though a slight preference for uP{x) > (i"(a;) and d^{x) < u^{x) at high 
X was noted. Changes in u"^ of up to ±10% were permitted. From conservation of quark 
number the percentage violation for d" is half that for u^. For both valence and sea 
quarks the consequent change in the proton distributions in the isospin- violating sets of 
partons is always 2% or less. We observed that the possible isospin violation was certainly 
sufficient to account for the NuTeV sin^ 9w anomaly, and that the slight preference in 
the type of violation for valence quarks is indeed in the correct direction to reduce this 
anomaly. We also discussed the possibility that s ^ s, noting that other studies have 
indicated a small inequality, the most recent of these again being in the correct direction 
to help resolve the NuTeV sin^ 9w anomaly. 

We note that the most significant theoretical errors come from possible corrections to the 
fixed-order DGLAP framework. Hence within the region of applicability, the conservative 
partons are the most reliable, and the variation of partons, and of physical observables, under 
changes in Xcnt a-nd Qcuf some indication of the theoretical uncertainties, both inside and 
outside the conservative domain. The most reliable extractions of as{M^) therefore follow from 
the conservative fits and are equal to 

a^^"(M|) = 0.1165 ± 0.002(expt) ±0.003(theory), (37) 

^NNLO^^I^ = 0.1153 ±0.002(expt)±0.003(theory). (38) 

In the NLO case this is considerably below our default determination, showing that the increase 
in as is mimicking theoretical corrections beyond NLO DGLAP. 

As additional tests of theoretical stability, we examined W and Higgs cross sections at the 
Tevatron and the LHC At the Tevatron we only sample partons in the conservative domain. 
Nevertheless the NLO predictions vary significantly as Xcnt and Q^^^ are changed due to the 
re-adjustment of the partons. The NNLO predictions are much more stable, implying less 
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theoretical uncertainty at NNLO. Indeed, the estimated total theoretical and experimental 
uncertainty of about ±2% on aw at the Tevatron offers an attractive and precise luminosity 
monitor. 

For non-central W and Higgs production at the LHC, we probe partons below the conser- 
vative Xcut — 0.005. For an, which samples the gluon not too much lower than x — 0.005, 
we still have reasonable stability, especially at NNLO. On the other hand aw samples quarks 
further below x = 0.005 and the NLO prediction can vary by up to 20% with different choices 
of cuts. This implies that the theoretical uncertainty at small x at NLO is rather large. How- 
ever, at NNLO, aw is much more stable (varying by about 3%), suggesting that the theoretical 
uncertainty has been considerably reduced by the inclusion of NNLO splitting and coefficient 
functions. We hope to confirm that this is still the case when the complete NNLO splitting 
functions become available. Assuming that this is so, the total theoretical and experimental un- 
certainty at the LHC is about ±4%. Therefore again it can serve as a good luminosity monitor. 
Work on resummations may be able to reduce the theoretical uncertainty still further. 

Therefore we conclude that theoretical uncertainties can be dominant in some kinematic 
regions, especially when the physical quantity probes partons at small x and/or small Q^. 
There seems to be considerable advantage in working at NNLO, as compared to NLO, but for 
real precision a few observables may have to await theoretical developments in low x and/or 
low physics. 
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